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ABSTRACT 

Method A method for rapidly forming a stochastic model of Gaussian or related type, 
representative of a porous heterogeneous medium such as an underground reservoir, 
10 constrained by data characteristic of the displacement of fluids. The method 

—It-comprises construction of a chain of realizations representative of a stochastic model 
(Y) by gradually combining an initial realization of (Y) and one or more other 
realization(s) of (Y) referred to as a composite realization, and minimizing an objective 
function (J) measuring the difference between a set of non-linear data deduced from 

15 said -the combination by means of a simulator simulating the flow in the medixmi and 
said — the data measured in the medium, by adjustment of the coefficients of the 
combination. The composite realization results from the projection of the direction of 
descent of the objective function, calculated by the flow simulator for the initial 
realization, in the vector subspace generated by P realizations of (Y), randomly drawn 

20 and independent of one another, and of the initial realization. During the-optimization 
stag e, the chain thus obtained is explored so as to identify a realization that allows te 
minimiz e m inimizing t he objective function (J). In order to sufficiently reduce this-the 
objective function, sequentially constructed chains constructed one after the other a re 
explored by taking as the initial realization the optimum realization determined for the 

25 previous chain. The method may be used for 

ApplicationG : notably development of oil reservoirs for oxamplo . 
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BACKGROUND OF THE INVENTION 
FIELD OF THE IN\^ENTION F ield ot the Invention 

The present invention relates to a method for more rapidly forming more rapidly a 
stochastic numerical model of Gaussian or a_related type, representative of a porous 
5 heterogeneous medium (such as a hydrocarbon reservoir for example) calibrated in 
relation to data referred to as dynamic data , characteristic of the displacement of fluids 
in the medium such as, for example, production data (pressures obtained from well 
tests, flow rates, etc.). 

The method according to th e invention finds applications in the sph e r e of 
10 und e rground zon e modelling intend e d to g e n e rat e r e pr e s e ntations showing how a 
c e rtain physical quantity is distributed in a zone of the subsoil (permeability notably), 
be s t compatible with observed or m e asur e d data, in order for e xampl e to favour the 
developm e nt th e r e of. 

BACKGROUND OF THE I>fVENTION Description of the Prior Art 

15 Optimization in a stochastic context consists in determining realizations of a 

stochastic model which meet a set of data observed in the field. In reservoir 
engineering, the realizations to be identified correspond to representations, in the 
reservoir field, of the distribution of carrying properties such as the permeability or 
porosity. These realizations form numerical reservoir models. The available data are, for 

20 example, isolated permeability or porosity measurements, a spatial variability model 
determined according to isolated measurements or data directly related to the fluid 
flows in an underground reservoir, ire rthat is pressures, breakthrough times, flow rates, 
etc. The latte¥ -data are often not linearly related to the physical properties to be 
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modelled. A randomly drawn realization is generally not in accordance with the whole 
of the data collected. Coherence in relation to the data is integrated in the model by 
means of an inverse procedure rSee 

- Tarantola, A., "Inverse probl e m theory Problem Theory - Methods for data fitting 
5 Data Fitting and mod e l paramet e r estimation M odel Parameter Estimation *\ Elsevier 
Science Publishers, 1987. 

The simplest technique is therefore the trial-and-error method. This approach 
conGists in r andomly taking takes realizations until a realization meeting the data is 
obtained. It affords the advantage of conservation of the spatial variability of the model 
10 but it-requires a-prohibitive calculating calculation t ime. It is therefore very rarely used 
in practice. 

An option that is often preferred is based on the-gradients calculation. The gradients 
methods allow to modify m odifying a n initial realization in a direction of search which 
is estimated from the gradients. The modifications are applied iteratively until data 
15 calibration is considered to be acceptable. The gradients methods are attractive because 
of their efficiency. However, they are no longer suitable when the number of 
parameters, ^re rthat is the number of permeability and porosity values forming the 
numerical model, is large. Besides, they do not allow to modify modifying the 
realizations while respecting the spatial structure of the stochastic model. 

20 More recently, a geostatistical parameterization technique has been introduced to 

constrain, by gradual deformation, the stochastic realizations to data on which they 
depend non-linearly. ft -This technique is described in French p atents FR:-2,780,798 and 
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FR.2,795,841 filed by o^the 



issignee. and in the following publications. 



notably : 



- Hu, L.Y., 2000, Gradual doformati eft -Deformation and itorativo calibration Interative 
Calibration of Gaussian rclatod GtochaGtic models Related Stochastic Models : Math. 

5 Geology, Vol.32, No. 1 , 

- Le Ravalec, M. Et al., 2000, The FFT moving av e rag e Moving Average (FFT-MA) 



gonorating — Generating and conditioning — Conditioning G aussian — simulations 

Simulations : Math. Geology, Vol.32, No.6, 

10 - Hu, L.Y., Blanc, G. And Noetinger, B. (2001) : Gradual d e formati en -Deformation and 
iterativ e calibration Interative Calibration of sequential stochastic Gimulations Sequential 
Stochastic Simulations . Math. Geology, Vol.33, No.4. 

This method has been successfully applied in various cases, notably from data 
obtained from oil development fields, as described in the following documents : 

15 - Roggero, F. et al., 1998, Gradual d e formati en — Deformation of continuous 
g e ostatistical mod e ls Continuous Geostatistical Models for history matchin gH istorv 
Matching , paper SPE 49004 : Proc. SPE Annual Technical Conference and Exhibition, 
New Orleans, 

- Hu, L.Y. et al., 1998, Constraining a roocrvoir facioa model Reservoir Facies Model to 
20 dynamic data u s ing Dynamic Data Using a gradual d e formation m e thod Gradual 
Deformation Method, paper B-01 : Proc. 6th European Conference on Mathematics of 
Oil Recovery (ECMOR VI), 8-1 1 September 1998, Peebles, Scotland. 




enerator : An efficient num e rcial m e thod Efficient Numercial Method for 
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As remindod d escribed in detail hereafter, the gradual deformation method allows te 
gradually modifV gradual niodification a realization of a stochastic model of Gaussian 
or a_Gaussian-related type while respecting the spatial structure of the model. 

Stochastic optimization 

5 Let = /2°'%. ., /m°'0 be the field data and / = (/n /z,..., /m) the 

corresponding responses simulated for a realization z of a given stochastic model Z. In 
general, the responses / = (/,, /a,..., /m) are obtained by solving iiumerically the direct 
problem. Thus, if z represents a permeability field, data / can be pressure 
measurements. In this case, they are simulated from a flow simulator. The goal of a 
10 stochastic optimization is to produce realizations of Z which reduce the differences 
between the observed data and the numerically simulated corresponding responses. 
These differences are measured by the following objective function : 

Coefficients are weights assigned to data which are functions of 

15 realization z. In this sense, minimization of the objective function is a problem with 
several variables. 

Let N be the number of grid cells or of components forming realization z. N is often 
very large (10^*- 10^. It is therefore very difficult to perform an optimization directly in 
relation to the components of z. Furthermore, realization z, even modified, must remain 
20 a realization of Z. The gradual deformation method allows these difficulties to be 
overcome. 

Random search from the gradual deformation method 
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W e now consid e r a A random field Z is now considered that can be transformed 
into a Gaussian random field Y. The gradual deformation technique allows to oonotruct 
construction of a continuous chain of realizations by combining an initial realization yo 
of Y with another realization ui, referred to as complementary, of Y, ui being 
5 independent of yo (Figure la). The combination coefficients are for example cos(t) and 
sin(t), and the combined realization meets the relation : 

y(t) = yo cos t + ui sin t 
where t is the deformation parameter. 

Another realization chain construction technique consists in combining combines 
10 the initial realization no longer with one, but with P complementary realizations Up 
(p=l,P) of Y (Figure lb). The coefficients of the combination are such that the sum of 
their squares is 1 . 

Once the chain is formed, it can be explored by varying the deformation parameter t 
and on e tries while trying to identify, from among all the realizations of this chain, the 

15 realization which minimizes the objective function. This minimization is performed in 
relation to t. Parameterization according to the gradual deformation method allows te 
reduc e r eduction of t he number of dimensions of the problem from N to 1 , where N is 
the number of values that constitute the field to be constrained. Furthermore, the sum of 
the combination coefficients squared being 1 , the optimized realization y(topt) still is a 

20 realization of Y : it- The optimized relation follows the same spatial variability model as 
all the realizations of Y. 

However, if the exploration of the r e alizations realization space is restricted to a 
single chain, euf -the p ossibilities of sufficiently reducing the objective function are 
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greatly limited. The above procedure therefore has to be repeated, but with new 
realization chains. These realization chains are constructed successively by combining 
an initial realization which is here the optimimi realization determined at the previous 
iteration, with a complementary realization of Y, randomly dmwn each time. Thus, at 
5 iteration /, the continuous realization chain is written as follows : 

yi (0 ^ yi-i cos t + ui sin t . 

yi.] is the optimum realization defined at iteration I A and the Ui are independent 
realizations of Y. The latter are also independent of yo. This formulation implies that the 
realization chain corresponds to a h3T)erellipse of dimension N. 

10 Minimizing the objective function in relation to t allows to improve improving or at 

least to preserve calibration of the data each time a new realization chain is explored. 
This iterative minimum search procedure is continued as long as data calibration is not 
satisfactory. The random side of the method lies in the fact that, upon each iteration, the 
complementary realization is drawn at random. In fact, the direction of search that is 

15 followed from the optimized realization at the previous stage is random. The direction 
of search, for a given chain and firom the optimum realization defined above, is : 



This direction of search only depends on ul. Furthermore, ul being independent of 
the complementary realizations already generated ui, U2,..., um and also of yo, the 
20 direction of search at the start of each new chain is orthogonal to the tangent defined for 
the previous chain at the same point (Figure 2). Although it may seem appropriate to 
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select a direction of search orthogonal to this tangent, there is an infinity infinite number 
of possible orthogonal directions. 

Experience shows that, after several iterations, the new directions of search no 
longer contribute significantly to the decrease of the objective function (Figure 6), 

5 It has also been considered combining to combine the initial realization not only 

with one, but with several complementary realizations. In this case, the number of 
deformation parameters increases rand it is equal to the number of complementary 
realizations involved in a gradual combination. Although the optimization process is 
then more flexible, several parameters have to be managed, which is not easy. Besides, 

10 such a process is not necessarily more efficient because it can depend on the execution 
of a larger number of direct flow simulations. 
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SUMMARY OF THE INVENTION 

The mefhod according to the invention allows more rapid formation of a numerical 
model representative of the distribution of a physical quantity in a porous heterogeneous 
medium such as an underground zone (oil reservoirs, aquifers, etc.), constrained by data 
5 collected in the medium (dynamic data characteristic of the displacement of fluids in the 
medium), collected by measurements (in production, injection or observation wells for 
example) or previous observations. 

The method according to the invention has applications for underground zone 
modelling that generates representations showing how a physical quantity is distributed 
10 in a subsoil zone (permeability notably), best compatible with observed or measured 
data, in order for example to favor the development thereof. 

It comprises an iterative process of gradual deformation wherein an initial 
realization of at least part of the selected model of the medium is linearly combined 
with at least a second realization independent of the initial realization, the coefficients 
15 of this linear combination being such that the sum of their squares is 1, and an objective 
function measuring the difference between a set of non-linear data deduced from said 
the combination by means of a medium simulator and gaid -the g eologic and dynamic 
data is minimized by adjusting the coefficients of the combination, the iterative process 
being repeated until an optimum realization of the stochastic model is obtained. 

20 ¥he -With the m ethod is characteriz e d in that the rate of gradual deformation to the 

optimum model representative of the medium is accelerated by selecting as the second 
realization to be combined with the initial realization at least one composite realization 
obtained by selecting beforehand a direction of descent defined as a function of the 
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gradients of the objective function in relation to all the components of said -the initial 
realization. 

The composite realization is obtained for example by linear combination of a set of 
independent realizations of the model, the coefficients of the combination being 
calculated so that the direction of descent from the initial realization y is as close as 
possible to the one- realization defined by the gradients of the objective function in 
relation to all the components of the initial realization. 

Optimization is for example carried out from a deformation parameter which 
controls the combination between the initial realization and the composite realization. 

In cases where said -the combination affects only part of the initial realization, the 
iterative process of gradual deformation is applied to a-Gaussian white noise used to 
generate a Gaussian realization and the derivatives of the objective function with 
respect to the components of the Gaussian white noise are determined. 

According to an implementation mode, the initial realization is combined with a 
certain number M of composite realizations, all obtained by composition from Pm 
independent realizations of Y, the optimization involving M parameters. 

In other words, the method oGsontially comprises a new gradual combination 
scheme which takes account of accounts for the information provided by the gradients 
at the initial point of any chain of realizations. Construction of a chain always starts 
from an initial realization and from a set of complementary realizations, all independent 
and coming from the same stochastic model. The initial realization is however not 
directly combined with the complementary realizations. The complementary 
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realizations make it possible to explore the realizations space in different directions. 
These directions are not equivalent : some allow to get nearer better approach to the 
optimum. At this stage, a realization referred to as composite realization is elaborated 
by combining the complementary realizations only. A chain of realizations is then 
5 created from the initial realization and from this composite realization. This chain, like 
the chain proposed in the basic case of gradual deformation, can be explored by means 
of a single deformation parameter. 

The composite realization is constructed go as t o propose provide a direction of 
search along the chain as close as possible to the direction of descent given by the 
10 gradients. As mentioned above, all complementary realizations are not equivalent : the 
composite realization takes the best of each complementary realization. 

The method allows to r e ach r eaching the formation of the numerical model 
representative of the medium more rapidly. 

BRIEF DESCRIPTION OF THE FIGURES 

15 Other features and advantages of the method according to the invention will be 

clear from reading the description hereafter of an application given by way of non 
limitative example, with reference to the accoinpanying drawings wherein : 

- Figures la, lb show gradual deformation schemes (referred to as GDM) that are 
already known, 

20 - Figure Ic shows the gradual deformation scheme (GBC) corresponding to the method 
according to the invention. 
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- Figure 2 shows realization chains in a Euclidean space with N dimensions, where the 
tangent at the level of the optimized realization for a realization chain /-I (RCy.i) is 
orthogonal to the direction of search for the initial realization of the next realization 
chain / (RCi), 

5 ' Figure 3 shows the projection of the direction of descent v in the subspace U defined 
by the orthonormal base formed by P independent realizations (ui,..., Up), 

- Figures 4A to 4C show permeability distribution examples respectively for the 
reference realisation, the initial realization and the realization constrained to the 
pressure data, 

10 - Figures 5 A to 5E respectively show the variations as a function of time of the 
bottomhole pressures respectively simulated for the five wells (BHP-OBSl, BHP- 
OBS2, BHP-PROl, BHP-OBS3, BHP-OBS4) of Figures 4, respectively for the 
reference (data), initial and constrained (match) permeability distributions, and 

- Figure 6 shows the evolution (OF) of the objective function as a function of the 
15 number k of flow simulations performed, GDMl corresponding to an optimization 

carried out by combining an initial realization and a single complementary realization, 
GDMGBC3, to an optimization carried out by combining the initial realization and^a 
composite realization constructed from three complementary realizations, GDMGBCIO, 
to the optimization carried out by combining the initial realization and a composite 
20 realization constructed from ten complementary realizations, and GDMGBC30, to an 
optimization carried out by combining the initial realization and a composite realization 
constructed from thirty complementary realizations. 

DETAILED DESCRIPTIO N OF THE INVENTION 
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The method according to the invention allows to orient upon orienting each iteration 
of the construction of the realization chain in order to reach a desirable direction of 
search. The technique selected takes advantage of the information provided by the 
direction of descent defined by the gradients of the objective function J. it-The 
5 technique can be implemented by means of a numerical simulator of a type known to 
the man skilled in the art, such as the ATHOS or ECLIPSE simulators. 

Search oriented by the direction of descent 

W e considor at At this stage consideration is given to t he direction of descent 
(evolution towards a minimum value, whether local or not) defined by the gradients of 
10 the objective function J in relation to all the components of realization z. These 
gradients are deduced from sensitivity coefficients 9/m / 9zn : 

The problem of sensitivity coefficients calculation has been widely dealt with in the 
scientific literature. The following document can for example be referred to : 

15 - Sun, N.Z., Inverse problems Problems in groundwater niodellin g Groundwater 
Modelling. Kluwer Acad. Publ., Dordrecht, The Netherlands, 1994. 

In particular, the adjoint state technique allows to caloulato calculation of all of 
these-coefficients by solving the direct problem and its -the adjoint problem, within a 
time interval equivalent to twice the time required for solution of the direct problem. 

20 Since the goal is to integrate the information provided by the gradients in the 

gradual deformation method, w e hav e to go going b ack to the Gaussian realization is 
necessary which results from the transformation of z. Two cases can be distinguished 
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according to whether isolated conditioning data, ire rthat is measurements of z at certain 
points, are available or not. 

When no data are available, if y is the Gaussian realization obtained by 
transforming z, then : 

dJ dJ dz„ 



5 

In the opposite case, we assum e the assumption is t hat z is known at certain points 
a. Thus, the Gaussian realization s, deduced from the transformation of z, is a 
conditional realization of Y. To obtain s, we-genefate-a non-conditional realization y of 
Y is generated and condition conditioned i t-to the values known at a by kriging. The 
10 conditional realization is deduced from : 

s = s* + (y - y*). 

s* and y* are respectively calculated by kriging from the real data arid the data of y 
generated at the level of points a. We-lLcan then shew- be shown that : 

dJ dz„ 



dJ 



-,Vn ^ a 
dz ds 

0, \/n = a 



15 For continuous physical properties, dzn/dyn or dzn/dsn are calculated from the 

anamorphosis fiinction allowing to transform realization z into a Gaussian realization. 
When the physical properties considered are category or discrete properties, these 
derivatives do not exist. The gradients techniques then cannot be applied. 

These various relations are of direct interest if a realization just has to be deformed 
20 globally. On the other hand, if it is desired to deform k -the realization locally, the 
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gradual deformation method has to be appHed to the Gaussian white noise x used to 
generate y. In this case, an additonal stage is required : calculation of the derivatives of 
the objective function with respect to the components of the Gaussian white noise. 

To illustrate the calculation of these derivatives, wo suggest to concontrato 
concentration is made on the particular case where the non-conditional Gaussian 
realization y is obtained from the FFT-MA generator described in the article published 
by Le Ravalec et al. 2000 mentioned above. 

The basic principle of this FFT-MA (FFT-Moving Average) generator is to 
transform a Gaussian white noise x into a Gaussian realization y correlated from a 
convolution product : 

y = g*x. 

Function g results from the expansion of covariance function C such that C = g * g\ 
where g^ is the transpose of g. The derivatives of the objective function with respect to 
the components of the Gaussian white noise are : 

dJ _^ dJ dy^ 

The discrete expression for the convolution product leads to * , which 

implies dyi/dxn = gi.n- If wo introduc e t his formulation is introduced in the derivatives 
of the objective function, w e show it can be shown t hat : 

This formula expresses the fact that the derivative of the objective function with 
respect to the n* component of the Gaussian white noise is given by the n**^ component 
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of the field obtained by convoluting all the derivatives of the objective function with 
respect to the components of the Gaussian realization with the kernel of the covariance 
function. From the framework established for the FFT-MA generator, these derivatives 
are determined as follows. 

5 1 -Calculation of the Fourier transform of dJ I 9y, these derivatives being are 

obtained by means of the direct numerical simulator ; 

2- Multiplication of this Fourier transform with that of g that is provided by FFT- 
MA ; 

3- Calculation of the inverse Fourier transform of the previous product. 

10 The time required for calculation of these derivatives is negligible : it represents an 

additional time of 2/3 in relation to the simulation of a Gaussian realization by FFT- 
MA. 

Whatever the realization generator and the direct numerical simulator, we assume it 
is assumed hereafter that we can d e fin e a direction of descent from the derivatives of the 
15 objective functio n can be defined . If optimization of the objective function is performed 
in relation to this direction of descent only, the coherence of the realization in relation to 
the spatial variability model is generally destroyed. In the section hereunder, we 
integrat e t he information provided by the derivatives of the objective function is 
integrated in the gradual combination scheme. 

20 Taking account of the derivatives of the objective function in the gradual 

deformation process 



W e consid e r th e The r ealization chain yl(t) constructed from yO and from another 
realization u of Y (Figure la) is considered . Now, instead of using a complementary 
realization u as it is, we randomly draw P complementary realizations Up (p=l,2,...,P) of 
Y are randomly drawn and we substitute u is substituted for a combination of the Up=i,p 
5 (Figure Ic). This combination is not any combination : it follows the following 
construction mode : 

p p 

The resulting realization u is a realization of Y and it is independent of yO. 
According to Eq.l, u is also tiie direction of search calculated for chain yi(t) at the 
10 starting point yo. W e th e r e fore suggest to construct The method constructs u so that tfie 
direction u is as close as possible to the direction of descent given by the gradients at yo. 

We first d e fine the The space V of the vectors is first defined with N dimensions 
provided with the scalar product : 

15 yin and yjn are respectively the n*^ components of vectors y; and yj . Let U be-is 

defined as a subspace of V defined by the orthonormal base (ui,U2,.- jUp). The direction 
of search in U which is the closest to the direction of descent v is given by the 
projection of v in U (Figure 3) : 

p 

20 By normalizing this vector, wo obtain the desired direction u is obtained . The 

combination coefficients X of Eq.(2) are thus : 
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The realization u thus defined is referred to as composite realization. 

We hav e considered so So far the construction of realization chains by combination 
of the initial realization with a composite realizatio n has been considered . This 
technique can however be generalized to the construction of chains involving a certain 
number M of composite realizations, all obtained by composition from Pm realizations 

ZM p 
^ , . realizations being independent. Optimization then involves M 

parameters. This technique increases the degree of freedom in the optimization process, 

but i^requires management of M optimization parameters. 

Numerical example 

We construct a A svnthetic reservoir model is constmcted on which the method 
according to the invention is tested. 

The synthetic reference reservoir is shown in Figure 4A. It is a monolayer reservoir 
comprising 51x51 grid cells which are^lO m thick and 40 m in side. The permeability 
distribution is lognormal with an average of 200 mD and a standard deviation of 100 
mD. The logarithm of the permeability field is characterized by an isotropic spherical 
variogram aiid a correlation length of 480 m. The other petrophysical properties are 
constant : the porosity is 25 %, the total compressibility 5.10"^ bar'^ and the fluid 
viscosity 1 cP. A production well BHP-PROl, with a radius of 7.85 cm, free of any skin 
effect, is at the c e ntr e center of the reservoir ^4tand is surrounded by four observation 
wells (BHP-OBSl, BHP-OBS2, BHP-OBS3, BHP-OBS4) (Figures 4). A numerical 
flow simulation allows to obtain obtaining, for this reservoir, a set of reference data 
comprising the bottomhole pressures of the five wells (Figure 5). 
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The object of the inverse problem is to determine a reservoir model coherent with 
the pressure data, the permeabilities distribution being assumed to be unknown. Four 
optiinization processes are therefore launched, starting from the same initial realization 
(Figure 4B). For each process, w e consider a single optimization parameter is 

5 considered , ire rthat is the deformation parameter. The first process (GDMl) takes up the 
conventional gradual deformation scheme with construction of a realization chain using 
the initial realization and a single complementary realization. The other three processes 
(GDMGBC3, GDMGBCIO and GDMGBC30) illustrate the application of the method 
according to the invention : the chaiiis are in this case elaborated with the initial 

10 realization and a composite realization obtained from the combination of 3, 10 and 30 
complementary realizations. The composite realization is constructed as explained 
above (Eq.2) by integrating the information provided by the gradients of the objective 
function with respect to the Gaussian white noise. For each process, the evolution of the 
objective fimction (OF) as a function of the number k of flow simulations performed is 

15 shown in Figure 6. 

It can be observed that using the gradients and increasing the number of 
complementary realizations allows the objective function to decrease more rapidly for 
the same number of flow simulations. 



